Ap. J., in press (1996) 



Structure of Stationary Photodissociation Fronts 

B.T. Draine 

Princeton University Observatory, Peyton Hall, Princeton, NJ 08544, USA; 
■ draine@astro.princeton.edu 

o\ . 

CT\ , and 

H ! Frank Bertoldi 

Max-Planck-Institut fiir Extraterrestrische Physik, D-85748 Garching, Germany; 



0O 



9 L,: 
6 

N 



C3 



fkb@mpe-garching . mpg . de 
ABSTRACT 



> 
(N 
CO . 

, The structure of stationary photodissociation fronts is revisited. H2 self-shielding is 

discussed, including the effects of line overlap. We find that line overlap is important 
for iV(H2) ^ 10 20 cm -2 , with a factor-of-two suppression of pumping rates at column 

' densities iV(H2) ~ 3 x 10 20 cnr 2 . We compute multiline UV pumping models, and 



compare these with simple analytic approximations for the effects of self-shielding. 

The overall fluorescent efficiency of the photodissociation front is obtained for 
different ratios of x/ n u (where x characterizes the intensity of the illuminating 
ultraviolet radiation) and different dust extinction laws. The dust optical depth r p dr 
to the point where 50% of the H is molecular is found to be a simple function of 
a dimensionless quantity (fio depending on x/ n H, the rate coefficient R(T) for H2 
formation on grains, and the UV dust opacity. The fluorescent efficiency of the PDR 
also depends primarily on (j)Q for x ^ 3000 and bh <^ 10 4 cm~ 3 ; for stronger radiation 
fields and higher densities radiative and collisional depopulation of vibrationally-excited 
levels interferes with the radiative cascade. We show that the emission spectrum from 
the PDR is essentially independent of the color temperature T co i or of the illuminating 
radiation for 10 4 K <^ T co i or , but shows some sensitivity to the rotation- vibration 
distribution of newly-formed H 2 . The l-0S(l)/2-lS(l) and 2-lS(l)/6-4Q(l) intensity 
ratios, the ortho/para ratio, and the rotational temperature in the v=l and v=2 levels 
are computed as functions of the temperature and density, for different values of x/ n ~R- 
We apply our models to the reflection nebula NGC 2023. Apparent inconsistencies 
between published K-band and far-red spectroscopy of this object are discussed; we 
adjust the two sets of observations for consistency. We are able to approximately 
reproduce the (adjusted) observations with models having x = 5000, kh = 10 5 cm~ 3 , 
and a reasonable viewing angle. Further observations of NGC 2023 will be valuable to 
clarify the uncertain spatial structure of the emission. 
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Subject headings: Infrared: ISM: Lines and Bands - ISM: Reflection Nebulae - 
Molecular Processes - Radiative Transfer - Ultraviolet: ISM 

1. Introduction 

It is by now widely recognized that reprocessing of OB starlight in photodissociation fronts 
plays an important role in the overall energetics of star-forming galaxies, particularly extreme 
"starburst" galaxies. An important mechanism in this reprocessing is the absorption of ultraviolet 
photons by molecular hydrogen, followed by ultraviolet and infrared fluorescence, and, about 
15% of the time, by dissociation. The vibrational fluorescence process, first noted by Gould & 
Harwit (1963), has been investigated by a number of authors (e.g.: Black & Dalgarno 1976; Shull 
1978; Black, Porter & Dalgarno 1981; van Dishoeck & Black 1986; Black & van Dishoeck 1987; 
Sternberg 1986, 1988; Sternberg & Dalgarno 1989). 

A number of theoretical investigations have advanced our understanding of photodissociation 
regions, or "PDRs".[] Solomon (1965: see Field, Somerville, & Dressier 1966) apparently was the 
first to propose that destruction of interstellar H2 might be dominated by dissociating transitions 
to the vibrational continuum following permitted transitions to the 1 BT,+ state; the resulting 
photodissociation rate was first estimated by Stecher & Williams (1967), who called attention to 
the potential importance of H2 self-shielding. The H2 self-shielding process has been studied by 
Shull (1978), Federman, Glassgold & Kwan (1979), and de Jong, Dalgarno, & Boland (1980). 

Tielens & Hollenbach (1985a) discussed the overall thermal and chemical structure of 
photodissociation fronts, and applied their theoretical model to explain observed properties of the 
photodissociation region in Orion (Tielens & Hollenbach 1985b). Black & van Dishoeck (1987) 
examined in detail the fluorescent excitation of H2, and the resulting infrared line emission. 
Recently Abgrall et al. (1992) presented new data relating to the photodissociation of H2, and 
carried out a more accurate treatment of the self-shielding process in order to reexamine the 
structure of the H/H2 transition zone. Recent detailed self-shielding calculations for a cloud of 
moderate density subject to strong UV illumination (Le Bourlot et al. 1992) and for diffuse cloud 
conditions (Heck et al. 1992) show that in some regions line overlap significantly affected the H2 
photodissociation rate. 

In a previous paper (Bertoldi & Draine 1996) we discussed the structure of coupled ionization- 
dissociation fronts, and concluded that they were in many cases expected to be propagating at 
significant velocities, calling into question the interpretation of some regions, including Orion, 



1 We use "PDR" for "photodissociation region". Note, however, that some authors use these initials for "photon- 
dominated region". 
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which have been based on models of stationary photodissociation fronts. In the course of our 
investigation we have reexamined the self-shielding of H2 . 

In this paper we revisit the problem of H2 self-shielding. We confirm that line overlap can 
often be important, and we identify the region in parameter space where this occurs. We develop 
an approximate method to allow for the effects of line overlap in a statistical fashion. We apply 
this method to compute models for photodissociation fronts including the effects of line overlap; 
we also improve somewhat upon previous treatments by applying the results of recent studies 
of the wavelength dependence of dust extinction. We find a simple analytic description of the 
H2 self-shielding function which appears to give a good approximation to the results of detailed 
self-shielding calculations, including line overlap. 

We construct models of stationary photodissociation fronts for different densities, 
temperatures, and intensities of incident FUV radiation. We examine how observable properties, 
such as the 1-0S(1) surface brightness, and the l-0S(l)/2-lS(l) and 2-lS(l)/6-4Q(l) line ratios, 
depend on the model parameters. Two indicators of the excitation mechanism (fluorescent vs. 
shock) and the gas density are the ratio of the ortho- (odd J) to para-H2 (even J) level populations 
and the rotational temperature within a vibrational level. We compute these quantities as a 
function of gas density and FUV illumination. 

The famous reflection nebula NGC 2023 is a well-studied example of fluorescent emission by 
H2. We attempt to reconcile published observations of this object. The available observations of 
NGC 2023 appear to be consistent with a photodissociation front with wh ~ 10 5 cm~ 3 irradiated 
by FUV radiation with intensity (relative to the Habing field) x ~ 5000. 



2. Computation of Level Populations 

2.1. Equations of Statistical Equilibrium 

Let the N vibration-rotation levels of the ground electronic state of H2 be designated by the 
index I = 1, N. Then 



dni 
dt 



= Rn u n(R)5i + ^ 



Aim + Pirn + ncC l 



tin 



n, 



A m i + n c C, 



ml 



^(/W + E^L,/) , (1) 



rfn(H) 
dt 



-2Rn H n(K) + 2J2 



Pdiss,m + E nc ^diss 



(2) 



The notation is as follows: R is the rate coefficient for formation of H2 via grain catalysis; Si is 
the fraction of H2 formed on grain surfaces which leaves the grain in level I; A\ m is the Einstein 
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A coefficient for spontaneous decay m — ► i (Ai m = for i£ TO < E{); fli m is the effective rate for 

(c) 

transitions m — ► / via ultraviolet pumping; C}' is the rate coefficient for transitions m — > i due to 
collisions with collision partner c G {H, H2, He, H + , e}; Cj^ is the rate coefficient for collisional 
dissociation out of level m by collisional partner c; fldiss,m is the rate of photodissociation out of 
level m. 

The effective transition rate from level m to level I via electronically excited states u due to 
ultraviolet pumping is 



u -ft-tot,u 

where Q u i is the rate for photoexcitation to electronically-excited level u out of level /, and the 
total rate of spontaneous radiative decay of level u is 



i 

here A vc ^ u is the rate of spontaneous radiative decay from level u to the vibrational continuum of 
the ground electronic state, resulting in dissociation. The dissociation probability for level u is 




(3) 




(4) 



Pdiss,u — 




(5) 



and the rate for photodissociation from level m is 




(6) 



u 



The dissociation probability averaged over photoexcitations out of level m = X(v,J) is 




(7) 



a 



It is convenient to define the diagonal elements An, Cu, and flu to be 




(8) 




(9) 




(10) 



Eq. (|T|) then becomes 
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where 

Dl m = A lm + Am + E "c^S > ( 13 ) 
c 

Ddiss,m = 0diss,m + nc ^diss,m 1 (^) 
c 

We include the 299 bound states of H2 with J < 29 in our calculations. 

2.2. Radiative Rates 

We designate the vibration-rotation levels of the ground electronic state X l T,+ by X(v, J), 
and the first 3 electronically-excited states B 1 ^, C 1 !^ and C 1 !!" by B(v,J), C + (v,J), and 
C-(v,J). 

Energy levels E u , transition probabilities A[ u , and dissociation probabilities Pdiss,u have 
been published by Abgrall & Roueff (1989) and Abgrall et al. (1992,1993a,b). We used data 
generously provided by Roueff (1992), covering levels up to J = 29. As discussed by Abgrall et 
al. (1992), these new data include dissociation probabilities Pdiss,u for the C~^(v, e/) levels which 
are considerably larger than the dissociation probabilities of Stephens & Dalgarno (1972), which 
were used for both C + (v,J) and C~(v,J) in earlier detailed modelling (e.g., van Dishoeck & 
Black 1986; Black & van Dishoeck 1987). As a result, we find somewhat larger overall dissociation 
probabilities. 

The transition probabilities A[ u and dissociation probabilities Pdiss,u are assumed to be 
independent of the energy of the photon responsible for photoexcitation to level u: the branching 
ratios for radiative decay out of level u are assumed to not depend on whether the photoexcitation 
to level u was due to a photon which was precisely resonant (A = A; M ) or far out on the damping 
wings of the transition (e.g., A = 1.001A/ U ). While this approximation is usual, it should be 
remembered that it is probably not exact. 



2.3. Collisional Rates 

We include inelastic collisions with H, He, H2, H + , and e~ , with collisional rate coefficients 
as described in Draine & Bertoldi (1996). For H-H2 collisions, inelastic cross sections have been 
computed by Mandy & Martin (1993,1996) using semiclassical trajectory calculations; Martin 
& Mandy (1995) have provided fits to the resulting rate coefficients for 450 < T < 20000 K. In 
the present application it is necessary to extrapolate to temperatures T < 450 K. The Martin & 
Mandy fitting function is 

,2 /4500K A 
log 10 (o-v) =a + bz + cz - d I — 11 , (15) 
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where z = log 10 (T/4500 K). Even for downward transitions, the coefficient d is typically large 
enough that this function declines very rapidly with decreasing temperature at low temperatures. 
Furthermore, comparison of semiclassical trajectory calculations and quantal calculations indicate 
that the former approximation substantially underestimates the inelastic cross sections near the 
energy threshold (Lepp, Buch & Dalgarno 1995). Accordingly, we use eq. (|l5|) for T > 6 = 600 K, 
but for T < 9 we take 

l2 /9000K , 4500KT\ 
log w {av) = a + bz + cz - d I — 1 J . (16) 

This expression joins smoothly to eq. (15) at T = 6, but declines less dramatically for T < 
Since we have extrapolated in a rather arbitrary manner, it is clear that quantal calculations with 
an accurate potential surface are needed to provide better estimates of the low temperature H-H2 
rate coefficients. 

H + is treated as a species with a fixed abundance xh = n(H + )/nn = 1 X 10~ 4 ; we assume^] 
n e /nu = 3 x 10~ 4 . 



2.4. H2 Formation on Grains 

As discussed elsewhere (Draine & Bertoldi 1996), we assume that H atoms recombine on 
grain surfaces to form H2 at a rate per volume 



( dn(K 2 ) \ 
V dt Ji 



R n H n(H) , (17) 



form 

with 



R = 6x 10~ 18 T 1 / 2 cm 3 s" 1 ; (18) 

this adopted value of it! is 2/3 of the value adopted by Black and van Dishoeck (1987) for most of 
their model calculations. For T = 70 K this rate is ~ 1.7 times larger than the value of R inferred 
by Jura (1975) from Copernicus observations of H2 in diffuse clouds. 

The rovibrational distribution 5(v, J) of the newly-formed H2 is taken to be of the form 
(Draine & Bertoldi 1996) 

5(v, J) = 5(0, 0)g n (J)(l + v) exp[-E(v, J)/kT f ) , (19) 



2 Even so, the resulting rate coefficients for (0, 3) -> (0, 1) and (0,2)->(0,0) at T = 200 K are only 2.9 x 10~ 14 cm 3 s" 1 
and 1.4 x 10 -14 cm 3 s _1 , smaller than the best estimates of Lepp et al. by factors of 30 and 100, respectively. 

3 For a photoionized metal abundance n(M + ) /tih ~ 2 x 10~ 4 (primarily C + ), a hydrogen ionization rate 
Ch = 8 x 10~ 18 (T/10 2 K) _0 ' 64 (riH/40cm _3 ) s _1 (due to either cosmic rays or X-rays) would sustain a hydrogen 
ionization fraction xn = 1 x 10 -4 . 
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where g n = 1,3 for even, odd J. We note that this distribution ( |i~9|) does not include a rotational 
degeneracy factor (2 J + 1) and does include a factor (1 + v); these deviations from a thermal 
distribution function are intended to enhance the populations of high-t; states relative to high- J 
states. In the present paper we take T f = 5 x 10 4 K, in which case the newly-formed H2 has an 
ortho/para ratio = 2.78, (v) = 5.3, (J) = 8.7, (J 2 ) 1 / 2 = 10.7, and a mean vibration-rotation 
energy of 2.89 eV. 



3. Radiation Field 

The efficacy of H2 self-shielding depends upon the distribution of H2 over the different levels 
X(v, J), which in turn depends upon the density and temperature of the gas, and the incident 
radiation field. 

In neutral regions, pumping of H2 is primarily effected by far- ultraviolet photons in the 
1110 — 912A range.0 Habing (1968) estimated the intensity of interstellar starlight at A = 1000A 
to be \u\ = 4 x 10 -14 erg cm -3 . We will consider various radiation fields u\, and will characterize 
the intensity of each at 1000A, relative to Habing's estimate, by the dimensionless factor 

y = 1 Jl 000 A (20) 

X ~ 4 x 10- 14 ergcm- 3 1 ' 

In Table 1 we list spectral forms which have been considered by various workers. For each 
spectrum we list the photon flux in the 1110 — 912A interval, 

11 10 A 

F= [ h^XuxdX , (21) 



912A 

the logarithmic derivative d\i\.u v / dhw evaluated at 1000A, and the color temperature Tcolor: 

the 

temperature of a blackbody having the same logarithmic derivative at 1000A. 

Habing (1968) estimated the average radiation field at 1000, 1400, and 2200A. We note that 
Habing's values can be fitted by 

\u\ = (-25Alj/6 + 25A 2 ; /2 - 13A 3 /3) x 10~ 14 erg cm" 3 , (22) 



4 Photons shortward of 912A are of course absorbed by atomic hydrogen. H2 in the v = levels has only weak 
absorptions longward of lllOA: the longest wavelength absorption out of the J — and 1 levels is the Lyman 
0-0P(l) line at 1110. 066A, with an oscillator strength fi u — 0.00058 (as compared to the much stronger 7-0R(l) (i.e. 
v u = 7, Vi = 0, J u — Ji = 1, Ji = 1) line at 1013. 434A, with fi u — 0.020). While there are weak absorptions out 
of v = 0, J = 2,3,4,5, ... longward of 11 10A (Lyman 0-0P(2),P(3),P(4),P(5),... at 1112.5, 1115.9, 1120.3, 1125.5A..., 
with f lu = 0.00070, 0.00074, 0.00076, 0.00076..., it is fair to assume that the bulk of UV pumping of H 2 is due to 
photons in the 1110-912A interval. 
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where A3 = A/10 3 A; this is the spectrum listed as "Habing" in Table 1. 

Draine (1978) approximated UV starlight by the relatively soft radiation field 

Xu x = 4 x l(T 14 xA3 5 [31.016A! -49.913A 3 + 19.897] erg cm" 3 , (23) 

with x = 1-71. 

Roberge, Dalgarno, & Flannery (1981) used a radiation field similar to (^)- This radiation 
field was employed by van Dishoeck & Black (1986) in their modelling of diffuse clouds, and by 
Black & van Dishoeck (1987) in their models of photodissociation fronts. As seen in Table 1, this 



spectrum is somewhat harder than eq. (23) 



As an example of the radiation field expected in a photodissociation front near an OB star, 
we consider the power-law spectrum u v oc v~ 2 : 

vu v {x = 0) = Xu x (x = 0) = 4 x 10" 14 xA 3 erg cm" 3 . (24) 

This spectrum has T co i or = 29000 K, corresponding to the spectrum of a B0 star. We use this 
spectrum in most of the model calculations; as we show below (§|6.2|, and Fig. 14), the properties 



of the PDR are insensitive to the spectrum for color temperatures T co i or ;> 10 4 K, so that the 
models presented here have general applicability to PDRs illuminated by stars of spectral type AO 
and earlier .0 

Pumping rates and dissociation rates for optically-thin H2 in various rotation-vibration states 
are given in Table 2 for the radiation fields (p4|) and (]23|). Also given in Table 1 is the probability 
Pdiss{v,J) that photoexcitation out of level X(v,J) will be followed by dissociation, and the 
probability (p re t(v, J)) that photoexcitation out of X(v, J) will be followed by direct spontaneous 
decay back to the original level X(v,J). Since there are many possible UV transitions possible 
out of any given level X(v,J), both Pdi ss and {p re t) depend on the shape of the illuminating 
spectrum, but Table 2 shows that for unshielded H2 the changes in {p^ss) or Pret are typically 
only a few percent when the shape of the radiation field is changed. Note, however, that changes 
in the radiation field due to self-shielding can lead to larger changes in, for example, {pdiss)i as 
the relative importance of inherently weak absorption transitions increases. Indeed, for complete 
photodissociation fronts we find (see § |5.2|) that (pdiss) ~ 0.15 averaged over all UV absorptions. 

Shull (1978) has noted that when the UV intensity is high, photoexcitation out of 
vibrationally-excited levels can affect the "infrared cascade" following UV pumping. For example, 
with no shielding the rate of photoexcitation out of the X(l, 3) level is 3.30 x 10 _10 xs" 1 for the 
spectrum (p4|). Since the total rate for spontaneous decay out of X(l, 3) is 8.35 x 10~ 7 s _1 , we see 
that in the optically-thin part of the photodissociation front (A^H^) ^ 10 14 cm~ 2 ) UV pumping 
will significantly interfere with the "infrared cascade" for \ k5 2500. Our calculations explicitly 
include UV pumping out of all levels X(v, J). 



5 Late B-type stars of course do not have large luminosities in the 1110-912A interval, and can only produce large 
X values close to the star. 
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4. Shielding of H2 



4.1. Shielding by Dust 



The wavelength-dependence of interstellar extinction has been discussed by Cardelli, 
Clayton & Mathis (1989, hereafter CCM) and O'Donnell (1994), and found to be very nearly 
a one-parameter family, where the parameter may be taken to be Ry = Ay/E(B — V). The 
extinction at A > 9000A appears to be "universal" , with an extinction cross section per H nucleus 
cw(9000A) = 2.4x 1(T 22 cm 2 . The "average" extinction curve has N H /E(B-V) = 5.8x 10 21 cm" 2 
(Bohlin, Savage & Drake 1978) and Ry = A V /E(B - V) = 3.1; the CCM extinction curve 
then predicts an extinction cross section per H of a ex t (1000A) = 2.62 x 1(T 21 cm . In dense 
molecular regions, however, the extinction appears to typically have larger values of Ry = 4 — 6 
(Mathis 1990), with a less steep rise in the ultraviolet. We will consider two cases: Ry = 3.1 
(for diffuse clouds) and Ry = 5.5 (the value observed toward ^ 1 Ori C in the Trapezium, possibly 
representative of very dense clouds); for Ry = 5.5 CCM predict a ex t(W00A) = 8.16 x 10~ 22 cm 2 . 

We consider plane-parallel clouds, with the radiation incident normally on one cloud surface. 
We do not attempt to treat scattering in any detail, but simply assume that the effect of 
dust is to attenuate the incident radiation field by a factor exp(— t^x), where t^x = ^H^d,x\ 
Nu = ^V(H) + 2iV(H2) + iV(H + ) is the total hydrogen column density between the point of interest 
and the cloud surface, and a^x is the effective attenuation cross section at wavelength A. 

The dust albedo in the 1110 — 912A region is uncertain; a recent study of the reflection nebula 
NGC7023 near A = 1000A finds an albedo ~ 0.4, and a scattering asymmetry factor (cos 6) ~ 0.75 
(Witt et al. 1993). The effects of scattering are of course complicated and geometry-dependent, 
but if the albedo is 0.4 then 0.6 < a d /a ext < 1. At A = 1000A we take the effective attenuation 
cross section to be cr^iooo = 2 x 10~ 21 cm 2 for Ry = 3.1 and 6 x 10~ 22 cm 2 for Ry = 5.5 - about 
75% of the estimated extinction cross sections at this wavelength. At wavelengths other than 
1000A we take the attenuation cross section a^x = (^a/^-iooo^^iooO; where the ratio of extinction 
cross sections is given by the CCM extinction curve for Ry = 3.1 or Ry = 5.5. 



Each level X(y, J) has many allowed transitions to vibration-rotation levels of the B, C + , and 
C~ electronic states. Let subscripts I and u denote the lower and upper levels. Let a continuum 
radiation field be incident on the cloud, with specific energy density ux at the cloud surface. 
Neglecting line overlap, the rate of photoexcitation I — > u for an H2 molecule in level I is simply 



4.2. Self-Shielding Neglecting Line Overlap 




(25) 
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where TV; is the column density of H2 in level I, r^iy) is the "attenuation" optical depth due to 
dust, and W u i(Ni) is the dimensionless equivalent width for the transition / — > u: 

W ul (N t ) = / y [1 - exp(-A^)] , (26) 

where the photoabsorption cross section cr u i(X), given by the Voigt line profile function, depends 
upon the oscillator strength f\ u for the transition, the usual Doppler broadening parameter 
b = FWHM/^/41n2, and the intrinsic broadening Atot,u of the upper level. We use accurate 
numerical approximations for dW u i/dN[ and W u i(Ni) given by Rogers & Williams (1974). 



4.3. Self-Shielding Including Line Overlap 



For small column densities, equation (|25|) will give a good estimate of the effects of 
self-shielding. At large column densities, however, lines may be shielded in part by overlap with 
other lines, as noted by Black & Dalgarno (1977). Many previous treatments have neglected 
this overlap: Black & van Dishoeck (1987) cited Roberge (1981) and Sternberg (1986) as having 
shown line overlap to be negligible for normal dust-to-gas ratios. More recently, detailed transfer 
calculations (Abgrall et al. 1992; Le Bourlot et al. 1992) have found that line overlap can become 
important under some conditions. Here we reconsider the importance of line overlap. 

Exact treatment of the effects of line overlap requires explicit calculation of the full radiation 
field as a function of frequency at each point. However, a simple estimate for the effects of line 
overlap may be obtained by first noting that the strong H2 absorption lines originate from the 
v = level and will therefore be almost entirely confined to the wavelength range 1110-912A ( sec 
footnote H|). Therefore for lines with A u ; > llloA we will neglect line overlap, and use equation 
(|2l 



We note that the wavelength range 1110 - 912A includes H Ly/3,7,<5, ... at 1025.72, 972.54, 
949.75, ...A. Since the PDR may develop large column densities iV(H), these lines may suppress 
the pumping rates for nearby H2 lines. 

For the lines with 912 < X u i < lllOA, we note that the total equivalent width of the 
ensemble of lines, including H lines and the effects of line overlap, will asymptotically approach 
W max ~ ln(1110/912) ~ 0.2 at large column densities. If we now treat the lines "democratically" - 
equally subject to the effects of line overlap - then we may approximate the effects of line overlap 
by considering a modified equivalent width W u i defined by 

= exp (-^/ w — ) for 912 < Kl < 11101 ' ( 27 ) 

where 

15 

w = E E w "i + E w ( m - n ) for 912 < A «< < llloA ; ( 28 ) 



l u n=3 
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the prime indicates that we restrict the sum to transitions with 1110 < X u i < 912A. For n ^ 10 
the H Lyman damping wings are not important, and the Lyman series lines blend together to form 
a continuum for n £ 39(3 kms" 1 /b) 1 / 3 , or A £ (911.76 + 0.60(6/3 kms^ 1 ) 2 / 3 )!. We will somewhat 
arbitrarily consider the H Lyman lines individually up to n = 15; the HI — 15 line is at 915. 83A. 
The damping constants Atot for the np levels were obtained by summing over A coefficients for 
the different decay channels; the latter were computed using radial integrals tabulated by Green, 
Rush, & Chandler (1957). We will see below that absorption by the H Lyman lines has only a 
small effect on the H2 pumping rates. 

It is easily seen that 

15 

W = Y, + E ^( H1 " n ) = W ™ax [1 - eM-W/Wmax)} , (29) 

U l n=3 

so that W — ► W max as W — > do. To allow for line overlap, therefore, we take the pumping rates to 
be 

(ul(Ni) = exp(-W/W max ) exp [-r d (X ul )\ , (30) 

with W u i given by eq. @ W given by eq. (|2§), and W max = 0.2. 

The pumping rate and photodissociation rate out of level I are simply 

(pump,l ^ t (ul i (31) 
u 

(diss,l ^ ' CulPdiss,u 

(32) 

u 

The effective pumping rate and photodissociation rate per FLj molecule are just 

(pump = [™(H 2 )] 1 ^ (pump,l n l , (33) 
I 

(diss = KH 2 )] -1 ECdiss,Z rt i i ( 34 ) 
where n; is the number density of FLj in level I. 



4.4. Approximations for Self-Shielding 

For many purposes it is useful to have simple analytic approximations for the dependence 
of the pumping rate on N 2 = iV(H2) and dust extinction t^iooo = -^h^iooo- Since the dust 
extinction does not vary strongly over the 1110 - 912A range, the pumping rate is approximately 

(pump{N 2 ) « fsMeid(N 2 )e- T ^ ooo („(0) (35) 
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where f s hieid{N 2 ) represents the self-shielding effect of the H2. Here we put forward two 
approximations for f s hield(^2)] their accuracy will be assessed later by comparison with detailed 
calculations. 



A very simple power-law approximation is provided by 

for N 2 < 10 14 cm" 2 
for N 2 > 10 14 cm 




fshield{N 2 ) — { ^ Ar /lnl4 -2\-0.75 at ^ -tnlA „ m -2 

(36) 



The N 2 75 dependence is steeper than the N 2 °' 5 behavior expected for heavily saturated lines, 



but we shall see below in §5.2 that it approximates f shield quite well over a large range of column 
densities. Note, however, that this approximation behaves unphysically in the limit N 2 — * 00: 
if there were no dust, the total pumping rate per area J °° ( pU mpdN 2 — » 00. This unphysical 
behavior corresponds to a failure to require that the total equivalent width be limited to some 
value W max ~ 0.2. We can correct this deficiency with a more complicated expression: 

fshieldiNz) = n °' 96 / t V2 + n ;°^ 05 exp [-8.5 x HT 4 (1 + xf 5 } , (37) 
(1 + x/b^y (l + x) u -° L J 

where x = N 2 /5 x 10 14 cm -2 , and b§ = 6/10 cms - . This somewhat complicated functional 
form has been constructed so that in the absence of dust the total pumping rate per area 
/o°° CpumpdN 2 w \F f° r t ne u v oc v~ 2 spectrum (p4|). It is also designed to reproduce the complete 
transition from optically-thin to extreme saturation with line overlap. We shall see below that 
equation ([37]) is significantly more accurate than equation (|3^), yet is nevertheless straightforward 
to evaluate numerically and hence useful for numerical modelling. 

Since the pumping rate may be written Cpump = F u dW u /dN 2 = vF u dW/dN 2 , (where 
F v = cu v /hu) and (pump = ( P ump{0) f shield, we obtain 

W(N 2 ) = Ainu p dN' 2 f shl eid{N' 2 ) (38) 

r Jo 

= Ahii/ 0.947(l+ •° 117 ^ -exp[-8.5x ^((l + x) - 5 -!)]! , (39) 
I l + x/6 5 L V 7JJ 



thereby demonstrating that our assumed shielding function ( |37| ) corresponds to an equivalent 
width W ~ Alni/ ~ 0.2 in the limit N 2 — > 00. Neglecting variations in the dissociation probability 
{Pdiss), we may approximate the photodissociation rate by 

(diss ~ fshieid(N 2 )e- T ^ooQ iss (o) . (40) 



5. Structure of Stationary Photodissociation Fronts 
5.1. Model Assumptions 



We now consider models of stationary photodissociation fronts, in which the H2 level 
populations at each point are the steady-state level populations resulting from UV pumping, 
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Fig. 1. — Comparison of self-shielding calculated using the statistical treatment of line overlap 
(eq.^0]) and an "exact" treatment of line overlap, for T = 100 K LTE level populations. The H2 
is assumed to have a Doppler broadening parameter b = 3kms _1 . Models do not include dust 
shielding or absorption by CI. It is apparent that eq. (30) does an excellent job of accounting for 



line overlap. Also shown are self-shielding approximations of Federman et al. (1979) and Tielens 
& Hollenbach (1985a) (see text). The upper panel shows the dissociation probability (pdiss) as 
a function of position; the dotted line shows the dissociation probability averaged over all UV 
pumping events in the PDR. 

spontaneous radiative decay, collisional excitation and deexcitation, photodissociation, and H2 
formation on grains. 

In the models presented here we do not examine the thermal balance in the PDR. Theoretical 
estimates for the temperature structure of PDRs are very uncertain, primarily due to uncertainties 
in the grain photoelectric heating: compare the temperature profiles of Tielens & Hollenbach 
(1987), Burton, Hollenbach & Tielens (1990), and Bakes & Tielens (1994). We will simply assume 
the kinetic temperature in the PDR to vary as 

T = TT~ > ( 41 ) 

1 + T (2,1000 

where the initial temperature Tq will be treated as an adjustable parameter. This temperature 
profile has the expected property of declining at large optical depth, as the dust grains attenuate 
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Fig. 2.— Same as Figure 1, but for T = 300 K LTE populations.. 

the incident UV. The functional form is not intended to be realistic, but does allow us to explore 
the effects of varying the gas temperature in the region where H2 pumping takes place. 

We assume a plane-parallel geometry, with radiation propagating in the +x direction. The 
UV intensity at x = is specified by the parameter x- At each point we assume an equilibrium 
between H2 formation on grains and H2 photodissociation. The H2 level populations are similarly 
assumed to be in steady-state statistical equilibrium. If we define 

yi = 2n ; /n H (42) 

and require yi = 0, then, since ra(H) = nn(l — xh — J2h=i Vl)i the vector y\ must satisfy the system 
of N inhomogeneous equations 

N 

(2iW* - D lm )y rn = 2Rn u (l - x H )<5/ . (43) 

m=l 

As discussed in §2, we set xh = 1 x 10 -4 , and include the iV = 299 bound states of H2 with J < 29 
in our calculations. We solve for the level populations y m ,m = 1,...,N by LU decomposition, 
using routines LUDCMP and LUBKSB from Press et al. (1992). We will refer to these models as 
"UV-pumping" models, since the rovibrational distribution of H2 is strongly affected by the UV 
pumping, and, as a result: (1) there is a significant population in rotational levels J > 5 which 
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Fig. 3. — H2 self-shielding factor f shield as a function of 2V(H2), in a stationary photodissociation 
front with x/ n K = -01cm 3 and To = 200 K. The H2 is assumed to have a Doppler broadening 
parameter b = 3kms _1 . Dust absorption is included assuming diffuse cloud dust with Ry = 3.1. 
The solid line shows the dissociation rate computed from a multiline calculation with approximate 
treatment of line overlap, using eq. (^). The long broken line is the semi-analytic fit (^)- Two 
power-law fits are shown: iV^ ' 75 (eqj3§), and a square-root fit (see text). The upper panel shows 
the dissociation probability (paiss)] the dotted line shows the dissociation probability averaged over 
all UV pumping events in the PDR. 

would be negligibly populated in a LTE distribution with T exc ~ 100 K; (2) the rovibrational 
distribution of the H2 varies with depth in the cloud, since pumping effects are strongest near the 
cloud surface. 

For comparison, we will also consider models in which there is a balance between H2 
photodissociation and formation on grains, but where the distribution of the absorbing H2 over 
rotation-vibration levels is assumed to be a thermal distribution with a specified excitation 
temperature T exc = 100 K. These will be referred to as "LTE" models. 
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Fig. 4. — Same as Fig. 3, but for dense cloud dust with Ry = 5.5. 



5.2. Self-Shielding 

In order to test the accuracy of our statistical treatment of line overlap, we have carried 
out detailed ("exact") radiative transfer calculations to compare with the photoexcitation rates 
computed using our statistical treatment of line overlap [eq. (p0|)]. The "exact" calculations 
explicitly took into account the 900 strongest FUV absorption lines of H2, with a fine enough 
frequency sampling to resolve all important frequency structure in the radiation field at different 
depths. This was achieved through an adaptive frequency mesh of typically 10,000 grid points and 
spectral resolution up to 10 6 , similar to the resolution in the calculations of Abgrall et al. (1992). 
Another 3200 weaker lines were included using the equivalent width approximation atop of the 
spectrum of the stronger lines. Typically, these weak lines contribute less than 5% to the total 
pumping rate. The CPU-time-efficient transfer algorithm is part of a fully time-dependent PDR 
code that will be described in a forthcoming paper. 



For comparison with the approximate treatment introduced in §4.3, the H2 level populations 
were assumed to be given by LTE at excitation temperature T exc . Rather than examine the 
actual self-shielding function, we consider the suppression of the dissociation rate, Cdiss/Cdiss(0). 
We compare the dissociation rate obtained from our "exact" radiative transfer equation with Cdiss 
obtained from eq. ©, for T exc = 100 K (Figure 1) and T exc = 300 K (Figure 2). We see that 
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Fig. 5. — Same as Fig. 3, but for dense cloud dust and a stronger radiation field, with 
x/ n H = 0.1 cm 3 . 

the agreement is excellent at all column densities ./V2. Agreement is of course expected for low 
column densities N2 ^ 10 18 cm -2 where line overlap effects are negligible (so that treatment of the 
pumping as due to isolated lines is valid). This agreement persists, however, at column densities 
A2 > 10 20 cm -2 where line overlap effects cause f shield to decline more rapidly than N^ ' 5 ; the 
agreement is excellent even out to the largest column densities considered, N2 = 2 x 10 21 cm -2 , 
where line overlap effects suppress the pumping rates by a factor ~ 10. We have therefore 
validated the use of eq. ( |30| ) for subsequent calculations, in which we solve for the non-LTE H2 
level populations. 

For comparison with our multiline calculations, in Figures |l] and ^ we have plotted the 
self-shielding function for a single line in the approximation of Federman et al. (1979), with the 
two parameters r = 1.3 x 10 -3 and A^/td = 10 14 ' 6 cm~ 2 chosen to approximately match the 
total H2 self-shielding curve for T exc = 100 K. It is apparent that the self-shielding of a single 
line is quite different from that of the total ensemble of overlapping H2 lines. The effects of line 
overlap were estimated by de Jong, Dalgarno, & Boland (1980); their result was used by Tielens & 
Hollenbach (1985a) for N% ;> 2 x 10 cm . In Figures I and I we show the self-shielding function 
used by Tielens & Hollenbach. The de Jong, Dalgarno, & Boland formula evidently underestimates 
the effects of line overlap, by approximately a factor of ~ 2.4 at N{R2) = 10 21 cm -2 . Also shown 
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Fig. 6. — Total dimensionless equivalent width W contributed by lines in the 1110 — 912A interval, 
as a function of the H2 column density, for x/ n u = 0.01cm 3 , b = 3kms _1 , and dense cloud dust. 
Separate contributions of Lyman band, Werner band, and HI Lyman series (Ly/3 through Lyman- 
15) are shown. Line overlap is expected to suppress pumping rates by a factor ~ 2 when W(total) 
equals 0.21n(2) = 0.14; this occurs at iV(H 2 ) = 3 x 10 20 cm- 2 . 

in Figures [l] and || is the dissociation probability (pdiss)- While {pdiss} varies from 0.10 to 0.21, it 
is generally close to the average value of ~ 0.15 (dotted line). 

Figures || and ||] show the variation of the self-shielding factor for dissociation, (diss/(diss{0), 
with depth in a cloud with x/ n K = 0.01cm 3 , for two extreme assumptions concerning the dust: 
diffuse cloud dust with Ry = 3.1 (Figure ||), and dense cloud dust with Ry = 5.5 (Figure Q). 
Figure [5] shows f shield for dense cloud dust, but with an enhanced radiation field x/ n n = 0-1 cm 3 . 
We see in Figures |^-[| that at H2 column density JV2 = 1 x 10 21 cm -2 , the H2 self-shielding 
suppresses the photodissociation rate by a factor f shield ~ 6 x 10~ 6 . It is also interesting to 
observe in Figs. |3|-|5] the decline in (pdiss) at large values of N2- This decrease (not seen in the 
dustless models of Figs. |l| and |2|) is due to attenuation of the FUV by the dust, resulting in 
"reddening" of the radiation to which the H2 is exposed, decreasing the relative importance of the 
photoexcitations to higher levels u, with larger values of Pdiss,u- 

Also plotted on Figures ||-|5] is the "square-root approximation" /shield ~ 
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- f dlss (0) = 4.42xlO- 10 s- 1 for UV pumping, x=10 
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Fig. 7. — Self-shielding for 3 photo-dissociation fronts: T = 100 K LTE level populations (solid line), 
and UV pumping with \ = 1 (dotted) and % = 10 (broken) (with nn = 10 cm" " 3 and T = 10 2 K). 
It is seen that the three self-shielding factors are essentially identical for iV(H2) ^2x 10 17 cm" 2 . 

7 x 10~ 5 (10 20 cm~ 2 /iV2) a5 ; it is evident that except over the limited range 10 17 ^ N 2 ^ 10 20 cm -2 
the H2 self-shielding function is quite different from the asymptotic N^ ' 5 behavior which is 
often assumed (e.g., Hill & Hollenbach 1978; Sternberg 1988; Goldschmidt &: Sternberg 1995; 
Hollenbach & Natta 1995). The differences arise from two effects. For 10 15 <^ -/V2 <^ 10 17 cm -2 the 
actual self-shielding function declines with increasing N2 much more rapidly than A^ 2 -0 ' 5 while 
many of the strongest absorbing transitions become optically thick and enter the "flat" portion of 
the cur ve-of- growth. The column density range of this steep drop is wider than that of a single 
line because initially less important transitions (i.e. those with smaller oscillator strength and/or 
smaller populations) are still optically thin and increase in relative importance when the stronger 
lines become saturated. Only at N2 > 10 17 cm~ 2 have all important pumping lines reached the 
square-root part of the curve-of-growth. For N2 ^ 10 20 cm -2 f shield falls off more rapidly than 
N 2 ~ ' 5 as the result of line overlap. The function 7 x 10" 5 (10 20 cm" 2 /^) 0,5 underestimates f shield 
by a factor ~ 20 at N2 = 10 15 cm" 2 , and overestimates f shield by a factor ~ 4 at ./V2 = 10 21 cm -2 . 

The importance of line overlap is shown by Figure ^, where we plot the total equivalent width 
W summed over Lyman and Werner series for H2, plus the H Lyman series from Lyman/? through 
HI — 15. Treating the frequencies of the lines in the 1110 — 912A band as independent random 
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Fig. 8. — H2 fractions in stationary, plane-parallel photodissociation fronts for bh = 10 2 cm -3 and 
To = 200 K, for selected values of x/ n n (cm 3 ) and dust with Ry = 3.1 {ad ,1000 = 2 x IO" 21 cm 2 ) 
and Ry = 5.5 (<7d,iooo = 6 x 10~ 22 cm 2 ). A > 912A radiation with u v oc v~ 2 is propagating in the 
+x direction at Nr = 0. A^h is the total column density of H nucleons. Self-shielding of the H2 is 
computed for 27983 lines using eq. ( |30"| ) with W max = 0.2. 

variables, our statistical treatment of line overlap, with J 'shield ^ ex P( — W/0.2), gives a factor of 
2 suppression of f shield when the total equivalent width W = 0.21n(2) = 0.14, which occurs for 
N 2 ss 3 x 10 20 cm" 2 (see Fig. |) 

If a power-law fit to the self-shielding factor /shield is desired, the 3 ^ 4 power law of equation 
(|36|), shown in Figures seems most suitable; it succeeds in reproducing f shield to within a 
factor of 2 for 10 14 < N 2 < 5 x 10 20 cm" 2 . 

If a more accurate approximate representation of f shield is desired, we favor equation (|37|), also 
shown in Figures Q-^. This functional form does an excellent job in reproducing the initial rapid 
decline in f shield for 10 14 < N 2 < 10 17 cm" 2 , the "square-root" behavior for 10 17 < N 2 < 10 20 cm -2 , 
and the rapid fall-off due to line overlap for N 2 > 10 20 cm -2 . 

As already seen from comparison of Figures |l| and ||, the detailed dependence of /shield on N 2 
depends upon the rovibrational distribution of the absorbing H2. In Figure we compare f shield 
computed for: (1) an LTE (T exc = 100 K) H2 level distribution (see Figure [l]); (2) a non-LTE 
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Fig. 9. — A^IFj) vs. Ah for photodissociation fronts with diffuse cloud dust. Models are labelled 
by x/ n H(cm 3 )- Calculations were done for «h = 10 2 cm -3 , To = 200 K. Line overlap contributes 
more than a factor of two to the self-shielding above the horizontal line. To the right of the vertical 
line dust shielding contributes more than a factor of 2. Points on each curve show locations where 
the H2 fraction is 0.1, 0.5, and 0.9. It is seen that for diffuse cloud dust properties dust shielding 
is important whenever line overlap occurs. 



rovibrational distribution due to UV pumping with \ = 1 an d collisional excitation/deexcitation 
with bh = 10 2 cm -3 and To = 200 K (see Figure and (3) same as (2), but with x = 10 (see 
Figure ||). In order to bypass any questions related to dust opacities, the plots in Figure are 
for cases with zero dust extinction. We see that for 10 14 <^ N2 <^ 10 18 cm -2 f shield does show a 
dependence on the details of the H2 rovibrational distribution. When there is higher excitation, 
the pumping occurs via a larger number of lines, so that the H2 remains optically-thin a bit longer, 
and self-shielding is less effective until one enters the heavily damped regime at N% ^ 10 18 cm -2 . 
Overall, however, the self-shielding function f s hieid{^2) is relatively insensitive to the effects of UV 
pumping with \ ^ 10 3 . 

Figure || shows profiles of the fractional abundance of H2 as a function of total column density 
A^h for five different values of x/ n n ( = 10 -3 , 0.01, 0.1, 1, and 10 cm 3 ) and two different dust 
opacity laws (Ry = 3.1 and 5.5). As the radiation intensity to density ratio x/ n H is raised, the 
column density of atomic, dusty gas increases. 
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Fig. 10. — Same as Fig. |9], but for dense cloud dust. It is seen that line overlap can be important 
even when dust absorption plays a minor role (e.g., x/ n R = -01cm 3 and N# = 5 x 10 20 cm~ 2 ). 
Calculations were done with n H = 10 2 cm" 3 , T = 200 K. 



In Figures ^ and |l0| we show the run of iV(H2) vs. Nn for models parameterized by x/ n H- 
Models assuming diffuse cloud dust properties are shown in Fig. ^; dense cloud dust is considered 



in Fig. 10. The region in which line overlap is important is shown in each figure, as is the 



region where attenuation of the UV by dust is significant. It is seen from Fig. 10 that line 
overlap plays an important role in self-shielding. For example, a dense cloud dust model with 
x/ n H = 0.01cm 3 has line overlap suppressing the dissociation rate by factors of 2-4 at column 
densities 5 <^ Nn ^ 12 x 10 20 cm~ 2 where dust extinction attenuates the radiation field by less 
than a factor of two. 



The 1-0S(1) line originates from H2(u = 1, J = 3). In Figure 11 we show the column density 
of this excited state, divided by Xi as a function of the total column density Nr, so that one can see 
how deep into the cloud one must go to account for most of the 1-0S(1) line emission. For each case 
we indicate the "median" for 1-0S(1) emission: the location where 50% of the 1-0S(1) emission 
occurs on either side. For x/ n H = 10~ 3 cm 3 , the 1-0S(1) median is at JVh = 1.5 x 10 20 cm~ 2 , but 
for x/ n H = 0.1cm 3 , for example, the 1-0S(1) median is at A^h = 8.5 x 10 20 cm~ 2 (for Ry = 5.5 
dust). 



10 17 10 1B 10 19 10 20 10 21 10 22 
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Fig. 11. — Column density of H-2(v = 1, J = 3) divided by \ in stationary, plane-parallel 
photodissociation fronts, labelled by x/ n H (cm 3 ), for dust with R v = 3.1 and 5.5 . Nr is the total 
column density of H nucleons. Self-shielding of the H2 is computed using eq. ( |30| ) with W max = 0.2. 
Calculations were done for kh = 10 2 cm -3 , To = 200 K. Results are valid for \ ^ 2000 and 
hh <; 10 5 cm" 3 (so that UV pumping and collisional deexcitation out of vibrationally-excited levels 
are slower than spontaneous decay). For each case a vertical line indicates the "median" location 
for 1-0S(1) emission (see text). 



5.3. Width of the PDR 



To assess the importance of shielding by dust, it is useful to define the dust optical depth of 
the PDR, T p< ir, which we take to be the dust optical depth at A = 1000A measured to the point 
where 50% of the local hydrogen is molecular: 2n(H2)/nn = 0.5. In Figure 12 we plot T p d r as a 
function of the parameter (fio defined in Paper I: 
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(44) 



where the rate coefficient R for H2 formation on grains is evaluated at temperature To. For R 
given by eq. (18) this becomes 



5.15 
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(45) 
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Fig. 12. — The dust optical depth r p d r at A = 1000A between the ionization front and the point 
where 2n(H2) = n(H), as a function of the dimensionless parameter <j>Q. Results are shown for two 
different extinction laws. It is seen that T p d r may be approximated by r p d r ~ ln(l + 5(/>o). Numerical 
results are for To = 200 K. 

We see that shielding by dust is important {j p d T > 1) for fa > 0.5, but relatively unimportant for 
(f>o < 0.1. A good fit to the numerical results is provided by the simple fitting function (cf. Fig. 

HD 

r pdr ~ ln(l + 50o ) • (46) 



5.4. Pumping efficiency 

Because of the increased importance of dust for large 4>q, the efficiency of conversion of 
incident UV into H2 fluorescence is reduced. One measure of this is the "pumping efficiency", the 
ratio 

_ $ CyumpdN 2 

tpump — , \1<) 

where F is the incident flux of 1110-912A photons (cf. eq. 21). Since there is only a minor 
contribution to pumping by A > lllOA photons, e pU mp is essentially the fraction of the incident 
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Fig. 13. — The UV pumping efficiency € pump , the fraction of incident 1110 — 912A photons which 
are absorbed by H2 rather than dust, and e^os^i), the number of 1-0S(1) photons emitted per 
incident 1100 — 912A photon. Numerical results are for To = 200 K. The estimate ( |48| ) for e pump 
is shown as a solid curve. The estimate fl5C)| ) for e 1 _ s , (i) is shown for dense cloud dust with 
R(T ) = 8.46 x 10 _17 cm 3 s~ 1 for kh — ► (broken curve) and for wh = 10 2 , 10 4 , and 10 5 cm~ 3 
(solid curves). 



1110 — 912A photons which are absorbed by H2, rather than dust, so that the H2 dissociation 
rate/area and (for «h ^ 10 4 cm -3 ) the fluorescent surface brightness of the cloud are proportional 

to £ p umpF C< £ p umpX- 

In Figure |^ we plot e pump versus the parameter 0o, for models with To = 200 K but a 
variety of densities and both diffuse cloud dust and dense cloud dust. We see that e pump 0.5 
for 4>q <; 0.1. Strictly speaking, e pump is a function of both x/nnR and cr^f] but we find that the 
pumping efficiency e pump can in practice be approximated as depending on the single parameter 



6 It is clear, for example, that for aa — > (in which case 4>o — ^ 0) we should have e pump » 1 since all 1110-912A 
photons must eventually be absorbed by H2 (neglecting the H Lyman lines and CI absorption). However, for values 
of ad and <j>o of practical interest, we find that e purnp « 0.5 for <f>o <C 1. 
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0o j with the fitting function 

0.5 

e P um P « YTfo ( } 

providing a reasonable fit to our results (cf. Fig. Also shown in Fig. [l3| is the "efficiency" 
e i-05(i) f° r emission of the strong 1-0S(1) line. We define this in terms of the surface brightness 
1(1 — > 0S(1), 9 = 0) measured normal to the PDR (and including the effects of internal dust; see 
eq. II): 

_ 4vrJ(1^0S(l),g = 0) 1 
ei-os(i) = r 7T ™ ■ (49) 

Provided that bh <^ 10 5 cm" 3 (so that collisional deexcitation does not compete with spontaneous 
decay out of (1,3)) we find that e^os^i) may be approximated by 

0.04 1 

ei -° S{1) ~ l + 20o [1 + (x/2000) - 5 ] ' (50) 

showing that for 0o ^ 0.1 and x <J 2000, e 1 _ 5(i) ~ 0.08e pump : ~8% of the UV pumping events 
result in fluorescent emission of a 1-0S(1) photon. Thus, for nn <J 10 5 cm -3 , 

/(1^OS(1),0 = O)« - ^-ergcm^s-^r- 1 . (51) 

V V ; ' ' (l + 20 o )[l + (x/2OOO) - 5 ] 6 V ; 



Unlike e pump , we see from Fig. |13| that e 1 _ s , (i) is not determined solely by 0o, but also shows 
sensitivity to the actual values of «h and x when rin ^ 10 5 cm" 3 or x ^ 10 3 . This is because at 
high densities or high intensities, (1) the rotational population is changed, affecting the probability 
of populating the v = 1, J = 3 level, and (2) collisional deexcitation (at bh ^ 10 5 cm -3 ) or 
radiative pumping (at x ^ 10 3 ) may act to depopulate the v = 1, J = 3 level. For nn 10 5 cm" 3 



we find that eq. fl50|) provides a good estimate for the intensity of the 1-0S(1) line. 

6. Fluorescent Emission Spectra 

6.1. Effects of Internal Extinction 

The true column density in level (v , J) along a line making an angle 9 with respect to the 
normal is simply 

If 00 

N(v,J) = j-r / dxn(v,J) . (52) 



cos(6>) Jo 

The emission in vibration-rotation lines of H2 is 



1 f°° Ai he 

I(u^l,9) = — / dx— !i n u —exp[-N }i (x)a{X u i)/cos(e)}, (53) 

cos(fl) Jo 4-7T X u i 

where 9 is the angle of inclination of the front as seen by the observer, and we have 
allowed for absorption of the fluorescent emission by dust within the photodissociation front 
(N n (x)=f x dx>n R (x>)). 
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When a surface brightness I bs(v, J — > v' , J') is observed in an H2 transition v, J — > ?/, J', one 
may directly compute the "apparent" column density N app (v, J) in the upper level from 

AT ( T\ 4TT\ u lI obs (v,J -» l/,J') 

J) = ^j^ajo • (54) 

Because of extinction, however, the true column density is larger. It is usual to infer an "observed" 
or "dereddened" column density 

N obs (v, J) = N app (v, J) 10° AA * , (55) 

where A\ is the estimated extinction, in magnitudes, at the wavelength A of the observed transition 
v,J—> v' J'. 



6.2. Sensitivity of the PDR to Incident Spectrum 

PDRs may be produced by radiation from stars ranging from late B-type to early O-type, and 
the properties of the PDR will in principle depend upon the color temperature of the illuminating 
radiation as well as on Xi measuring the intensity at IOOOA. Fig. 14 shows ei_os(i) an d selected 



line ratios for a series of PDR models where x an d «h are held constant but the illuminating 
radiation is taken to be a dilute blackbody with color temperature T co i or = 10000 K (AO spectrum) 
to T co i or = 40000 K (04 spectrum). The models have % = 10 3 an d «h = 10 4 cm 3 - for this density 
collisional deexcitation is relatively unimportant for vibrationally-excited levels, so that the 
emission from these models should show maximum sensitivity to the spectrum of the illuminating 
radiation. We see, however, that the emission properties of the PDR are essentially independent 
of T C oior for T co i or > 12000 K. This insensitivity allows us to limit our study to models computed 
for power-law spectra with a = —2 (corresponding to T co i r = 29000 K) and apply them to PDRs 
produced by a broad range of star types. 



6.3. Sensitivity to S(v, J) for Newly-Formed H2 



In Fig. 14 we also investigate the sensitivity of the emission spectrum to the assumed 
distribution function S(v, J) of newly-formed H2. The solid curves are for models with "hot" H2, 
with Tf = 5 x 10 4 K, and the broken curves are for models with "cold" H2, with Tf = 2000 K. We 
see that the intensities of lines out of levels with low J [e.g., S(l) lines out of J = 3] are essentially 
unaffected by the assumed properties of 5, even for high vibrational levels [e.g., 7-5S(l)]. Lines 
out of levels with higher J values show increasing sensitivity to the assumed 5: S(5) lines out of 
J = 7 levels vary by about 40% between models with hot or cold H2, and S(9) lines out of J = 11 
levels are about a factor of ~ 2 stronger in models with hot H2 versus cold H2. 

It is interesting to note that the sensitivity to the assumed S(v,J) is primarily determined 
by the J value of the upper state, and seems not to depend particularly on the value of v. the 
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Fig. 14. — Selected line intensity ratios for x = 10 3 , nH = 10 4 cm -3 , versus the color temperature 
Tcoior of the incident radiation field. Curve labelled ei_os , (i) is the efficiency for 1-0S(1) line emission 
[eq.(|49"|)]. Solid and broken curves are, respectively, for newly- formed H2 which is rovibrationally- 
"hot" (T f = 5 X 10 4 K) and "cold" (T f = 2000 K). The PDR is assumed to have T = 500 K, dust with 
Ry = 5.5, Od.1000 = 6 x 10~ 22 cm 2 , and H2 formation rate given by eq. (18). Line ratios are seen to 
be essentially independent of T co i or for T co i or ^ 12000 K. Intensities of lines out of levels (v,J) with 
J = 3 [S(l) lines] are essentially independent of the value of Tf. Intensities of lines out of levels with 
J = 7 [S(5) lines] are increased by ~ 50% when Tf is varied from 2000 K to 5 x 10 4 K. Intensities 
of lines out of levels with J = 11 [S(9) lines] are a factor of ~ 2 stronger when the newly- formed 
H 2 is "hot" rather than "cold" . 



10-4S(5) line and the 1-0S(5) line each increase by about 50% as the newly-formed H2 is changed 
from "cold" to "hot". Unfortunately, as we show below, the rotational distribution of the H2 
depends on x> ^H, and the gas temperature; unless these other properties of the PDR are known 
quite accurately, it does not seem likely that we will be able to use observations of PDR emission 
spectra to determine the distribution function 5(v, J) of newly-formed H2. 



6.4. 1 0S(l)/2 1S(1) line ratio 



The l-0S(l)/2-lS(l) line ratio is frequently used to characterize H2 line-emitting regions. In 
Fig. [D] we show this intensity ratio (for PDRs viewed face-on) as a function of gas density riH, 
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Fig. 15. — l-0S(l)/2-lS(l) intensity ratio vs. density nn- The dust is assumed to have Ry = 5.5, 
0cZ,iooo = 6 x lCP 22 cm 2 , and R(T) from eq. (|l8|). Curves are labelled by To (see eq. ^jj). For pure 
fluorescence (n H < 10 5 cm -3 and T < 500 K) we have l-0S(l)/2-lS(l)=1.9; collisional effects 
increase the ratio by collisional excitation of v = 1, and preferential deexcitation of u = 2. 



and a number of values of the temperature To [see eq. (pl~|)]. Results are shown for two different 
values of the ratio x/ n H- At low densities and temperatures the line ratio (for u u oc v~ 2 , or 
Tcoior ~ 3 x 10 4 K) is I(1-0S(1))/I(2-1S(1))=1.9, reflecting excitation by radiative pumping, with 
negligible collisional deexcitation (or radiative excitation) out of vibrationally-excited levels. As 
seen from Fig. 14, the line ratio is insensitive to the value of T co \ or . 



As the temperature and density are increased, two effects occur: collisional deexcitation and 
collisional excitation. Collisional deexcitation out of the v = 2 (and higher) levels acts to reduce 
the intensity of 2-lS(l) emission; while collisional deexcitation also depopulates the v = 1 levels, 
the rate coefficients for deexcitation out of v = 2 levels (by collisions with H; Martin & Mandy 
1995) are significantly larger than those for deexcitation out of v = 1, hence collisional deexcitation 
tends to raise the l-0S(l)/2-lS(l) line ratio. 

If the temperature is high enough, collisional excitation of the v = 1, J = 3 level can also 
become important. For x/ n H = 0.01cm 3 and nn <^ 10 4 cm -3 , collisional excitation contributes 
~ 5% of the 1-0S(1) emission for T = 700 K, ~ 30% of the 1-0S(1) emission for T = 800 K, and 
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Fig. 16. — As in Fig. 15, but for the 2-lS(l)/6-4Q(l) intensity ratio. The line ratio is affected by 
the rotational distribution of the H2 (which explains why the line ratio depends upon To even for 
riH = 10 2 cm" 3 , and upon the effects of collisional deexcitation of (v, J) = (2, 3) and (6, 1). 



dominates the 1-0S(1) emission for To ^ 850 K. 

As x/ n H increases, collisional excitation becomes relatively less important (at fixed To) for two 
reasons: First of all, our adopted temperature profile (eq. [fil]]) has the highest temperature in the 
region where the H2 abundance is low. In the optically-thin region, the H2 fraction is proportional 
to nn/x, so that increasing x/ n H decreases the amount of H2 present in the hottest part of the 
PDR (e.g., where 0.8Tq < T < To). Secondly, when x/ n u 1S increased, T p d r increases (see Fig. 



|i2|) . For our assumed temperature profile (41), this implies a reduction in T at the location where 
2n(H2) = n(H), thereby decreasing the rate of collisional deexcitation. Therefore, when x/ n n is 
increased, collisional excitation of 1-0S(1) does not become important until higher temperatures: 
for x/ n K = 0.1cm 3 we see that for To = 1000 K the collisional excitation contributes only ~ 10% 
of the 1-0S(1) emission, and only begins to dominate for Tq ;> 1200 K. 
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6.5. 2 lS(l)/6 4Q(1) line ratio 

The 6-4Q(l) line (A = 1.6015/im) is sufficiently strong for spectroscopic measurement (see, 
e.g., Luhman & Jaffe 1996), and the 2-lS(l)/6-4Q(l) line ratio is a useful indicator of conditions 
in the region where UV pumping is taking place. Whereas the 2— 1S(1)/1-0S(1) line ratio can 
be affected by collisional excitation of (1,3) {E/k = 6951 K), the 2-lS(l)/6-4Q(l) line ratio is 
much less sensitive to collisional excitation because the (2,3) level has E/k = 12550 K). The 
2-lS(l)/6-4Q(l) ratio is therefore a relatively direct probe of collisional deexcitation of (6,1) vs. 
(2,3). " 

The dependence of the 2-lS(l)/6-4Q(l) line ratio on the density tt-h and temperature profile 
parameter Tq is shown in Fig. [R| Because the 2— 1S(1) line originates from J = 3 while the 
6— 4Q(1) line originates from J = 1, the line ratio depends on the gas temperature even at low 
densities where collisional deexcitation of the levels is negligible: the line ratio increases with 
increasing To at low densities. When the gas density is high enough that collisional deexcitation 
begins to compete with radiative decay of the (v, J) = (6, 1) level, the line ratio increases, since the 
collisional rate coefficients (Martin & Mandy 1995) are such that the (6, 1) level is more strongly 
affected than the (2,3) level. At the highest density «h = 10 6 cm~ 3 , we obtain very high values 
of 2-lS(l)/6-4Q(l) when T Q = 200 K (see Fig. [jjj). At these high densities the relative strengths 
of 2-lS(l) and 6-4Q(l) are primarily determined by the relative rates of collisional deexcitation 
of (2, 3) and (6, 4) by atomic H. We remind the reader that we use an uncertain extrapolation of 



the results of Martin & Mandy (1995) to low temperatures (see eq. (16). For the temperatures 
T < 200 K in our To = 200 K models, these extrapolated rates are very uncertain; we have included 
these models mainly to show that extreme values of the 2-lS(l)/6— 4Q(1) line ratio are possible 
at high densities. Reliable calculations of line ratios for high density, low temperature models 
will not be possible until accurate inelastic cross sections are available for low energies, requiring 
quantal calculations and an accurate potential surface. 



6.6. Models 

A PDR has hundreds of fluorescent emission lines which may be observable, but space 
limitations preclude discussion of other than the 1-0S(1), 2-lS(l), and 6-4Q(l) lines discussed 
above. Complete H2 vibration-rotation emission spectra may be obtained via anonymous ftp for 
the models listed in Table ^. Readers interested in H2 spectra for other values of x> n H> Tq, ad, 
Tf, or b should contact the authors. 
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Fig. 17. — Same as Fig. 15, but showing ortho/para ratio parameter 7 [see eq. (p9[)] for H2 in levels 
v = 1, 2 < J < 7. An H + fraction xh = 10 -4 has been assumed. 



7. Rotational Populations 

The rotational populations within a given vibrational level can be approximated by a thermal 
distribution, 

N(v, J)/gj = aexp(-E(v, J)/kT r ) , (56) 

where 



gj = (2J + 1) for even J ("para-H 2 ") 
= 3(2J + 1) for odd J ("ortho-H 2 ") 



(57) 
(58) 



and T r is the apparent "rotational temperature". Because ortho- and para-H2 may not be fully 
equilibrated, we allow separate values of a, denoted a ort ho and a para , for the ortho- and para-H2 
states, but we fit a single rotational temperature T r for both ortho- and para-H2 within a single 
vibrational level. We define the ratio 



7 = 3a ortho /a 



para 



(59) 



to characterize the ortho-para ratio: from our definition of 7, we see that for a thermal distribution 
we would have 7 = 3, even at low temperatures. Note that 7 is not the actual ortho/para ratio 
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Fig. 18. — Same as Fig. [0], but showing rotational temperature of H2 in levels v = 1, 2 < J < 7. 

(i.e., the ratio of the total column densities of the ortho-H2 and para-H2), which goes to zero in 
thermal equilibrium as T —* 0. 

In Figure |l7] we show 7 determined by fitting eq. ( |56| ) to our model calculations; the fit is 
restricted to the v = 1 levels 2 < J < 7. For n H ^ 10 4 cm -3 and T 700 K, we find 7 2 ± 0.2, 
essentially independent of x/ n H- 

As discussed in §^~J|, the H2 is assumed to be formed on grains with an ortho/para ratio 
of 2.78, so an appreciably lower value of 7 for the v = 1 levels reflects ortho— >para conversion 
in the gas. We have assumed xh = n(H + )/nn = 10~ 4 in the present models. In regions with 
in > 10 -7 , the dominant ortho-para conversion process is reactive scattering with H + , with a rate 
coefficient ~ 2 x 10~ 10 cm 3 s _1 for H2(0, 1) + H + — > H2(0, 0) + H + . While this is much slower than 
photodissociation of H2 in unshielded regions, note that at the point where 2n(H2) = n(H) the H2 
photodissociation rate is (diss = Rn-u = 6 x 10~ 17 (T/100 K)°- 5 (?t,h/ cm _3 )s _1 , slow compared to 
the rate 2 x 10 _14 (?t.h/ cm -3 )(xH/10 -4 ) s _1 for (0,1)— >(0,0) conversion by H + . In regions where xh 
is very small, ortho-para conversion takes place via scattering by H atoms and off grain surfaces. 

As ran is increased above ~ 10 4 cm -3 , 7 increases and approaches ~ 3. Indeed, we see that 
this non-LTE environment can result in values of 7 > 3: for \ = 10 5 ) «h = 10 6 cm -3 , and 
Tq = 200 K, our best-fit 7 = 4.3. Note that this occurs under conditions where UV pumping can 
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Fig. 19.— Same as Fig. [18, but for levels v = 2, 2 < J < 7. 



compete with spontaneous decay from vibrationally-excited levels (cf. §3). The observed large 
value of 7 may result from the greater ability of ortho-H2 to self-shield. 

The rotational temperatures T r {v = 1) and T r {v = 2) of the v = 1 and 2 levels of H2 are 
shown in Figures ^ and 19, We see that optical pumping results in rotational temperatures 
1000 — 2000 K. There is a general trend for T r to increase with increasing x ( e -g-) compare lower 
and upper panels of each Figure), due to the increasing importance of optical pumping compared 
to radiative decay or collisional deexcitation. It is interesting to note that under some conditions 
an increase in To leads to a decrease in T r : for example, at nn = 10 4 cm~ 3 and % = 10 2 > changing 
To from 500 to 1000K causes T r to decrease from 1500 to 1000K. This is due to the large increase 
in H-H2 inelastic rate coefficients, making it possible for collisional processes to compete with 
radiative pumping, tending to bring T r toward the gas temperature T. 



8. NGC 2023 
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8.1. Observations 

The B1.5V star HD 37903, at an estimated distance D = 450 pc, is situated near the edge of 
the molecular cloud L1630, resulting in the prominent reflection nebula NGC 2023. NGC 2023 
displays strong H2 fluorescent emission, and its spectrum has been studied both in the infrared 
(Gatley et al. 1987; Hasegawa et al. 1987) and the far-red (Burton et al. 1992). Gatley et al. 
have mapped the 1-0S(1) emission, which peaks in an "emission ridge" SSE of the star. The 
brightest spot on the ridge is located about 78"S, 9"W of the star (Burton et al. 1989; Field et al. 
1994; Brand 1995), corresponding to a transverse distance of 5.3 x 10 17 (Z?/450pc) cm. Field et al. 
(1994) have obtained a high resolution image of NGC 2023 in the 1-0S(1) line; the brightest spot 
is 27"N, 22"E of the star, with a surface brightness slightly higher than on the southern emission 
ridge. The C recombination line emission peaks in the vicinity of the southern 1-0S(1) emission 
peak (Wyrowski & Walmsley 1996); C recombination line spectroscopy indicates a gas density 
rin ~ 10 5 cm -3 (Pankonin & Walmsley 1976, 1978). Since the only available far-red spectra of 
NGC 2023 were taken on the southern filament, our modelling objective will be to try to reproduce 
the observed H2 spectrum of the southern filament. 

A B1.5V star (T cff = 22000 K, L = 7600L©) radiates S uv ps 2.5 x lO^s" 1 in the 1110 - 912A 
range (Fitzpatrick 1995, based on Kurucz ATLAS9 model atmospheres); at a distance r from the 
star, 

_ S uv /4nr 2 ( 5xl0^cm\ 2 
X * 1.209 x MP cm-2 s -i = 6600 { r J ' (60) 

neglecting dust extinction. Harvey et al. conclude that most of the infrared-emitting dust is more 
than O.lpc from HD 37903. We therefore assume that there is only modest absorption by dust 
between the star and the H2 emission ridge, and seek to reproduce the observed H2 line spectra 
by a model of a plane-parallel stationary photodissociation region illuminated by radiation with 
X ~ 5000. The idealized geometry is shown in Figure |20|. 

The computed surface brightnesses will of course depend on the adopted value of the angle 
9. Since the highest fluorescent surface brightness will presumably be in a region where there 
is significant limb-brightening, we need to estimate the likely magnitude of this effect. A shell 
with inner and outer radii r, and r\ + Ar, containing material with emissivity j, would have 
a peak surface brightness (in the absence of internal extinction) j(8rjAr) 1 / 2 [l + Ar/^rj] 1 / 2 . 
Viewed face-on, a portion of the shell would have surface brightness jAr^, so that limb- 
brightening enhances the surface brightness over that of a plane-parallel slab viewed normally 
by a factor B = (8rj/Ar) 1 / 2 [l + Ar/2ri] 1 / 2 . For «h ~ 10 5 ns cm -3 and UV extinction cross 
section ciooo ~ 6 x 10~ 22 cm 2 , we would estimate the thickness of the fluorescing region to be 
Ar ~ (nnfTiooo) 1 ~ 1-7 x lO 16 ?!^ 1 cm. Taking ~ 5 x 10 17 cm we then estimate a peak geometric 
limb-brightening factor B ~ 15n^ 2 . We regard this as an upper limit to the value of (cos#) -1 . If 
the local radius of curvature of the PDR at the surface of a "filament" is smaller than the distance 
to the star, the peak limb-brightening factor will be smaller than this upper limit. 
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Fig. 20. — Slab geometry assumed for the H2 "emission ridge" observed for NGC 2023. 

A convenient way to compare with observations is to plot N b s (v, J)/gj vs. E(v,J), where 
E(v, J) is the energy and gj is the degeneracy of the (v, J) level. The "observed" N b s (v, J) values 
are obtained from eq. (|55|) J 7 ] We characterize the extinction by Ak, and have used the Cardelli, 
Clayton, & Mathis (1989) extinction curve for Ry = 5.5 to estimate A\/Ak- Burton (1993) has 
estimated the K-band extinction to be Ak = 0.3 mag; we assume a slightly smaller extinction 
Ak = 0.2 mag, which corresponds to extinction by dust with E(B — V) pa 1.37^4^ = 0.27 mag and 
Ay = 1.5 mag. The extinction at 0.9 and 0.7 /iin is estimated to be Aq$ = A.2Ak = 0.84 mag and 
A . 7 = Q.0A K = 1.2 mag. 

There are 3 levels (v = 3, J = 3, 4, 5) for which emission has been reported by both Hasegawa 
et al. (1987; hereafter H87) and Burton et al. (1992; hereafter B92). Unfortunately, after 
correcting for differential extinction, the column densities N b s (v, J) derived from the B92 far-red 
spectrum are a factor of ~ 10 — 20 larger than those derived from the K-band observations of 
H87! For example, the 3-2S(3) flux observed by H87, corrected for reddening by Ak = 0.2 mag, 
corresponds to iV b s (3, 5) = 2.5 x 10 cm" 2 . However, the 3-0S(3) flux reported by B92, corrected 



7 Note that the column densities N(v, J) (uncorrected for extinction) reported by Burton et al. (1992) are not 
consistent with their reported intensities, due to a multiplicative error (Burton 1995); the Burton et al. column 
densities must be reduced by a factor 0.36 to bring them into agreement with the reported intensities. 
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X = 5000, n H =10 5 cm- 3 T = 500K, cosfl = 0.20, A K =0 2 




E(v,J)/k (10 3 K) 



Fig. 21. — Column densities N(v, J) divided by the level degeneracy g(J) for the NGC 2023 
emission ridge. Filled symbols are derived from observations of Hasegawa et al. (1987) and Burton et 
al.(1992), corrected for extinction with Ak = 0.2 mag. Surface brightnesses reported by Hasegawa 
et al. (1987) have been multiplied by 6 to allow for beam dilution (see text); surface brightnesses of 
Burton et al. have been reduced by a factor 0.4 for consistency (see text). Open symbols are for a 
plane-parallel slab model with x = 5000, nn = 10 5 cm -3 , and To = 500 K, viewed at an angle with 
cos# = 0.2. Solid lines are least-squares LTE fits to the 2 < J < 7 populations in each vibrational 
level. For each vibrational level the best-fit values of the rotational temperature T r and ortho/para 
ratio parameter 7 are indicated. 

for reddening by ^4.7962^m = 5.08Ak = 1.02 mag, corresponds to N b s (3, 5) = 3.4 x 10 15 cm -2 - a 
factor of 13 larger than inferred from the K band spectra! Discrepancies of 21 and 15 are found 
for the (3,3) and (3,4) levels, although for these levels the reported uncertainties in the fluxes are 
greater. We will assume an overall discrepancy of 15 between the H87 and B92 column densities. 
The fluorescent emission is strongly peaked on the "emission ridge" (Burton et al. 1989, Brand 
1995), and the differences between H87 and B92 are in large part due to the fact that H87 used a 
diameter beam while B92 used a 27" x 1" aperture aligned along the emission ridge. We will 
assume that the H87 surface brightnesses need to be multiplied by a factor /3h87 = 6 to correct for 
beam dilution, but that the B92 fluxes should be multiplied by a factor /3b92 = /?H87/15 = 0.40 to 
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Fig. 22.— Same as Figure |l], but for model B, with X = 5000, n H = 10 5 cm" 3 , T = 900 K, and 
cos/9 = 0.2. 



correct for an apparent calibration error .^] 



8.2. Models for the NGC 2023 PDR 

We consider two stationary PDR models for the NGC 2023 emission ridge. Both have 
X = 5000, bh = 1 x 10 5 cm -3 , and a viewing angle 6 = cos -1 0.2 ~ 78°. We assume the temperature 
profile in the PDR to be given by eq. ( fH| ) and consider two values of To: model A has To = 500 K, 
and model B has Tq = 900 K. Both models have x/ n H = 0.05 cm 3 ; c/)q = 0.12 for model A and 



0.086 for model B. From Fig. 11, for both models we expect most of the vibrationally-excited 
H2 to be within the first Ah = 1.5 x 10 21 cm -2 , or 1.5 x 10 16 cm, corresponding to ~ 2" at the 
distance of NGC 2023. This is consistent with the observed narrowness of the "emission ridge" . 

The pumping efficiency e pump = 0.69 for model A and 0.74 for model B, somewhat larger 



Note that with the H87 intensities increased by a factor of 6, the B92 intensities would not require correction 
if there were no differential extinction between the K band and the far-red. This possibility does not seem likely, 
however, as even the dust within the PDR should contribute extinction with an effective Ak ~ 0.2. 
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than the estimates ~0.45 from eq. (f48|); the enhancement of e puTnp results from the substantial 
populations of H2 in excited rovibrational levels, so that UV absorption by H2 is spread over more 
lines, including many with A > IIIOA. The models succeed in reproducing most of the observed 
intensities to within a factor ~ 2. As we shall see, model B appears to be in fairly good agreement 
with the published observations for NGC 2023. 



8.3. 1 -0S(l)/2 -1S(1) Ratio 

Hasegawa et al. find l-0S(l)/2-lS(l) = 3.7 ± 0.6; from Figure 2 of Burton (1993) we infer 
a value of 3.0 for this ratio. These are both well above the calculated value for our model A, 
which has l-0S(l)/2-lS(l)=2.20, a value typical of UV pumping. Model B, with a higher gas 
temperature, has l-0S(l)/2-lS(l)=2.91; this is in part due to collisional excitation of the v = 1 
level, and in part to preferential collisional deexcitation of UV-pumped v = 2, as discussed in §|7[ 
Varying To between 800K and 1000K allows us to obtain l-0S(l)/2-lS(l) ratios between 2.54 and 
3.90. We consider Model B to be in good agreement with the observed l-0S(l)/2-lS(l) ratio. 



8.4. Rotational Temperatures and Ortho-Para Ratio 



In Figures Oil and 22 one sees that there is fairly good overall agreement between predicted 



and observed line intensities. Figures 21 and |22| show the least-squares fit to the rotational 
temperature T r characterizing rotational levels 2 < J < 7 for each vibrational level. H87 obtained 
a rotational temperature T r = 900 K for the v = 1, J = 2 — 4 levels. If we restrict our least-squares 
fit to these 3 levels, we obtain T r = 1217 K and 7 = 2.36 for Model A and T r = 1160 K, 7 = 2.46 
for Model B, in reasonable agreement with H87. 

H87 found T r ps 1500 ± 200 K for the v = 2, J = 2 - 5 levels. For these 4 levels, our 
least-squares fit gives T r = 1587 K for model A, and T r = 1650 K for model B. 

For the higher vibrational levels we also do a fairly good job of reproducing the relative 



populations of the different J levels, as can be seen for the v = 3, 4, 7, and 8 levels in Figures |21 
and |2^. 

Our ortho-para ratios are somewhat above the observed values. Burton (1993) finds that the 
observed line intensities are consistent with 7 = 2.1 ± 0.2 for all levels, whereas we find (for Model 
B) 7 = 2.46 for v = 1 (J = 2 - 4), and 7 = 2.46 for v = 2 ( J = 2 - 5). 
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I[l-0S(l)] = 3.29xl0- 4 erg crrr a s" 1 sr- 1 
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Fig. 23.— Emission spectrum for NGC 2023 model A (x = 5000, n H = 10 5 crcr 3 , T = 500 K) 
viewed from a direction with cos# = 0.2. The spectrum is convolved with a Gaussian response 
function with R = A/FWHMa = 2 x 10 3 . Extinction by dust within the PDR has been included, 
but no external dust has been allowed for. 



8.5. Spectra 

We have taken the predicted emission spectra for the above models and convolved them with 
a Gaussian with A/FWHM A = 2000 (FWHM„ = 150 km s" 1 ). The spectra include the effects of 
extinction within the emitting slab, but no additional foreground extinction is allowed for. In 



Figures |23| and 24 we show the predicted emission spectra for models A and B. The far-red lines 
derive from transitions from high-v states, for which the collisional deexcitation rate coefficients 
are much larger than for the low-f states. Note that the far-red emission is relatively stronger in 
model A than in model B - this is because there is more collisional deexcitation in model B due 
to the higher temperature. A list of all the H2 lines with intensities exceeding 0.1(1280 and 1174 
lines for models A and B, respectively) is available via anonymous ftp (see Table pi). 
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X=5xl0 3 , n H =10 5 cm- 3 , T =900K, cos0 = O.2 
I[l-OS(l)] = 4.04xlO" 4 erg cm" 2 s" 1 sr" 1 
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Fig. 24.— Same as Fig. |23|, but for model B with T = 900 K. 



8.6. Discussion 



Black &: van Dishoeck (1987) have previously put forward models for the fluorescent emission 
from NGC 2023. Their models had l\jy = 300 which (cf. Table 1) corresponds to x ~ 670 - 
a factor of 7.5 weaker than the x = 5000 favored here. The reported 1-0S(1) face-on surface 
brightness for their "reference model A" is 3.5 x 10 -5 erg cm -2 s _1 sr -1 , corresponding to an 
"efficiency" ei_os , (i) = -058. We have computed a model with x = 670, bh = 10 4 cm -3 , T = 80 K, 
and the same grain properties (extinction and H2 formation rate) as assumed by Black & van 
Dishoeck for their reference model A; we obtain an efficiency £i_o5(i) = 0.021, a factor of 2.7 
smaller than obtained by Black & van Dishoeck. Part of the discrepancy is attributable to the fact 
that Black Sz van Dishoeck did not allow for the effects of line overlap; if we modify our model 
to neglect line overlap (by taking W max — > 00 in eq. 50), the 1-0S(1) surface brightness increases 
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by a factor of 1.5. The remaining discrepancy, amounting to a factor of 2.7/1.5 = 1.8, may be 
due in part to the fact that we use improved (larger) dissociation probabilities for the C + levels, 
although this seems unlikely to explain the factor of 1.8. 

The Black & van Dishoeck models were developed with the aim of reproducing the surface 



brightnesses reported by H87. As discussed in §8.1, we now believe the peak 1-0S(1) surface 



brightness to be a factor ~ 6 higher than reported by H87, so the published Black & van Dishoeck 
models are not applicable to the bright bar in NGC 2023. 

Our modelling efforts lead us to conclude that if the dissociation front is approximately 
stationary, then the fluorescing H2 in the NGC 2023 emission ridge must have bh ~ 10 5 cm -3 , and 
appreciable limb-brightening (such as cos# = 0.2, as in our Models A and B), in order to reproduce 
the observed high surface brightness. Lowering nu (increasing x/ n u) would lead to somewhat 
lowered fluorescent efficiency, and by reducing the importance of collisional effects would make 
it difficult to reproduce the observed l-0S(l)/2-lS(l) line ratio, which requires a combination 
of collisional deexcitation from v = 2, and collisional excitation from v = to 1. On the other 
hand, increasing nu would result in increased collisional deexcitation of vibrationally-excited 
levels, making it difficult to reproduce the observed high surface brightnesses. It is reassuring 
that our density estimate nu ~ 10 5 cm~ 3 is in agreement with that deduced by Pankonin & 
Walmsley (1976, 1978) from C recombination lines, which should originate from the same gas. 
Jansen, van Dishoeck, & Black (1994) have used the observed HCO + 3-2/4-3 line ratio to infer 
nu ~ 5 x 10 4 cm -3 in this region. Fuente, Martin-Pintado, and Gaume (1995) have argued for 
nu « 10 6 cm -3 near the photodissociation front based on observations of CN, SO, and HCN 
rotational lines. 

Our best model (B) does not fully reproduce the observed line surface brightnesses (see Fig. 



22). There are a number of possible explanations: (1) The emission bar may include a range of 
densities, whereas in our modelling we seek to explain the observations with a single density. (2) 
The actual temperature profile may be rather different from that assumed here. (3) Line overlap 
effects will suppress particular photoexcitation transitions, whereas in our treatment we share this 
suppression among all of the transitions with 1110 < A < 912A. (4) As discussed in §^[^, the actual 
branching ratios following photoexcitation in the damping wings may differ from our adopted 
branching ratios, thereby altering relative line strengths. (5) Errors in our adopted collisional rate 
coefficients may be important, since the vibrationally-excited states may be collisionally-deexcited 
at the densities «h ~ 10 5 cm -3 assumed here. (6) The dissociation front in NGC 2023 may be 
propagating rapidly enough that the stationary approximation made here leads to errors. 

NGC 2023 shows structure on ~1" scales (Field et al. 1994). We need calibrated H2 
spectral-line images of NGC 2023 in a number of H2 fluorescence lines to test the accuracy of our 
models. 
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9. Summary 

The principal results of this paper are as follows: 

1. An approximate treatment of self-shielding [eq. (^0|)] is developed which treats line overlap in 
a statistical fashion. This treatment allows individual lines to be treated using the single-line 
equivalent-width approximation, but allows statistically for overall suppression of the 
continuum due to line overlap. Comparison with exact calculation shows this approximation 
to be accurate. 

2. Two simple approximations are provided for the self-shielding function for H2. The first, 



a simple power-law [eq. (36)], provides a fairly good fit to self-shielding over the range 



10 15 <N 2 < 10 21 cm . The second [eq. p7|)] is an analytic function which provides a very 
good approximation to the self-shielding including the effects of line overlap. These functions 
are recommended for use in future studies of photodissociation regions. 

3. The effects of line overlap become important for column densities N(H-2) ;> 10 20 cm -2 , 
suppressing the UV pumping rate by a factor of 2 at iV(H2) ~ 3 x 10 20 cm -2 (cf. Fig. ||). 

4. For dense cloud dust properties, the effects of line overlap become important while the dust 
is still optically thin for x/ n H ^ 0.05 (cf. Fig. [Tc|). 

5. The dust optical depth T pc i r at the point where 2n(H_2) = n(H) is primarily a function of a 
single parameter 0o oc (x/nnR)a %qq, defined in eq. (f44]) . The approximation (46) provides a 
good estimate for T pc i r . 

6. The "efficiency" of UV pumping of H2 in the photodissociation front, e purnp , is also a function 
of the single parameter <po, and is approximately given by eq. 



7. For rin <^ 10 5 cm 3 the efficiency ei_os(i) for 1-0S(1) emission is given as a function of cj>Q 
and x by eq. (jj(J). 

8. The emission from PDRs is quite insensitive to the color temperature of the illuminating 
radiation for color temperatures T co [ or ^ 10 4 K, or stars with spectral type AO or earlier. 

9. The emission spectrum from low J levels [e.g., 1-0S(1), or 3-2S(l)] is insensitive to the 
distribution function 8(v, J) of newly- formed H2, but the emission from higher J levels [e.g., 
1-0S(9)] is increased (by a factor ~ 2) when the newly-formed H2 is rotationally- "hot" . 

10. Observable properties of photodissociation fronts - including the l-0S(l)/2-lS(l) line ratio 



(Fig. 15), the 2-lS(l)/6-4Q(l) line ratio (Fig. 16), the ortho/para ratio (Fig. O), and 



the rotational temperatures (Fig. 18,19) - are computed for various values of x, nn, and 



temperature. Complete H2 vibration-rotation spectra are available via anonymous ftp. 
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11. The reflection nebula NGC 2023 is considered. We correct the H87 intensities for assumed 
beam dilution, and adjust the B92 intensities for consistency. The (adjusted) observations are 
approximately reproduced by a model with = 10 5 cm -3 and \ = 5000, with Tq ~ 900 K. 
The agreement is not perfect; possible explanations for the discrepancies are discussed. 

We are especially grateful to P.W.J.L. Brand and M. Walmsley for communicating recent 
observations of NGC 2023 in advance of publication, to E. L. Fitzpatrick for computing S uv for B 
stars, to R.H. Lupton for making available the SM plotting package, to P.G. Martin for making 
available H-H2 rate coefficients in advance of publication, and to E. Roueff for kindly providing 
us with molecular data for H 2 in computer-readable form. We wish also to thank M.G. Burton, 
E. van Dishoeck, D. Neufeld, M. Walmsley, D.T. Jaffe, and E. Roueff for helpful comments and 
discussions. 
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Table 1. Interstellar Ultraviolet Radiation Fields (A > 912A) 



Radiation Field 
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Table 2. Pumping and Dissociation Rates for Unshielded H2 



ri2 level 
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from eq. (|24|) . 
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Table 3. Library of Stationary PDR Models' 
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ouu 


1 

J. 


i ef; y i n -5 

J- .OO A J. U 


Hh3o 


10 4 


10 3 


6 


1000 


1 


2.20 x 10~ 5 


Lm3o 


10 5 


10 3 


6 


300 


1 


2.61 x 10~ 5 


Lw3o 


10 5 


10 3 


6 


500 


1 


2.69 x 10~ 5 


Mw3o 


10 5 


10 4 


6 


500 


1 


9.66 x 10~ 5 


Mh3o 


10 5 


10 4 


6 


1000 


1 


1.17 x 10~ 4 


Qm3o 


10 6 


10 4 


6 


300 


1 


1.32 x 10~ 4 


Qw3o 


10 6 


10 4 


6 


500 


1 


1.42 x 10~ 4 


Rw3o 


10 6 


10 5 


6 


500 


1 


3.86 x 10~ 4 


Rh3o 


10 6 


10 5 


6 


1000 


1 


4.66 x 10" 4 


n2023a 


10 5 


5000 


6 


500 


0.2 


3.29 x 10~ 4 


n2023b 


10 5 


5000 


6 


900 


0.2 


4.04 x 10~ 4 



a Models available via anonymous ftp from astro.princeton.edu, subdirectory draine/pdr, or 
at trttp : //www . astro . princeton . edu/ ^draine/ . All models have b = 3kms _1 , newly-formed 
H2 with Tf = 5 x 10 4 K, and incident UV with u v oc v~ 2 . 

b cf. eq.pl) 



c Parameter determining temperature profile via eq.(41). 
d 6 = viewing angle. 



